#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

sig.effects(measure = "FA", atlas = "glasser", dataset = "pnc")

sig.effects(measure = "FA", atlas = "glasser", dataset = "hcpd")

smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)

smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")

ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  scale_y_continuous(limits = c(-0.02, 0.0138)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)

#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
  
  #Read in participant-level data and model-level fitted values and derivatives 
  participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
  fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA

  #Connection spline plot with participant-level data
  trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
    geom_point(color = alpha(display_color, 0.5), size = .1) +
    geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
    geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted,  ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
    labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
    theme_classic() +
    scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    scale_y_continuous(breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  #Significant derivatives (rate of age-related change) plot 
  derivatives.plot <- ggplot(data = derivatives.df) + 
    geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
    scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color,  na.value = "white") +
    labs(x="\nAge", fill = "Derivative") + 
    theme_classic() +
    scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    theme(
        legend.position = "none",
        axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size = 6, family = "Arial", color = c("black"))) 
  
  allplots <- list(trajectory.plot, derivatives.plot)
  tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)

  return(tract.plot) 
}

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") 

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.06, 0.15)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") 

ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
  geom_point(color = c("#9c3a80"), size = 0.4) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
  scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Agematuration_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.94, height = 1.58)

cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)

cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late

matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey

ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Maturationalbrackets_PNC.pdf", device = "pdf", dpi = 500, width = 4.4, height = 4.4)

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))

devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.pnc

gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")

devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
  arrange(desc(term.correlation))

devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))

devpattern.neurosynth.hcpd

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)

ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 2.2, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  scale_y_discrete(
  labels = c("expectancy"=expression(bold(expectancy)),
  "monitoring"=expression(bold(monitoring)),
  "cognitive_control"=expression(bold(cognitive_control)),
  "response_inhibition"=expression(bold(response_inhibition)),
  "thought"=expression(bold(thought)),
  "recall"=expression(bold(recall)),
  "induction"=expression(bold(induction)),
  "object_recognition"=expression(bold(object_recognition)),
  "visual_perception"=expression(bold(visual_perception)),
  "gaze"=expression(bold(gaze)),
  "induction"=expression(bold(induction), parse=TRUE))) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_line(color = c("gray90")),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)

#Function to compute the number of overlapping terms when drawing term sets of size 20 at random for PNC and HCPD nulls
compute_term_overlap <- function(){
term.subset.pnc <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
term.subset.hcpd <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
num.overlap <- length(intersect(term.subset.pnc, term.subset.hcpd))
return(num.overlap)}

#Null distribution of overlapping term counts based on 10,000 random draws
set.seed(1)
overlap.number <- replicate(10000, compute_term_overlap()) %>% as.data.frame %>% set_names("overlaps")

sprintf("In %s out of 10,000 runs, the term overlap number was equal to or greater than 10", sum(overlap.number$overlaps >= 10))

sprintf("The permutation-based p value for these overlap counts is %s", sum(overlap.number$overlaps >= 10) / 10000)

summary(overlap.number$overlaps)

ggplot(overlap.number, aes(x = overlaps)) + 
  geom_histogram(bins = 12, fill = c( "#9c3a80"), color = c("white")) +
  theme_classic() +
  scale_x_continuous(limits = c(-1, 10), breaks = c(0, 2, 4, 6, 8, 10)) +
  labs(x = "\nNeurosynth term overlap", y = "Count\n") +
  geom_vline(xintercept = 10, linewidth = .1) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .1), 
  axis.ticks = element_line(linewidth = .1))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_overlap_null.pdf", device = "pdf", dpi = 500, width = 1.3, height = 1.1)

derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)

derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey

ggplot() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  17.5, na.value = "gray93", limits = c(16, 19), oob = squish) + 
  new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) 

S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)

ggplot() + 
  geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  180, na.value = "gray93") + 
   new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) 

cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_hcpd.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages 
colnames(temporalSAcorr_FA_glasser_pnc) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_pnc))) #set column names; 100 age x 10,000 posterior draw correlations df

deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws

cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")

ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
  geom_ribbon(fill = c(alpha("#323280", 0.4))) +
  geom_line(linewidth = .2, color = "black") +
  labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
  theme_classic() + 
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
  scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
  theme(
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.line = element_line(linewidth = .2), 
   axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivative-SAcorrelation_plot_pnc.pdf", device = "pdf", dpi = 500, width = 1.35, height = 1.03)

#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages 
colnames(temporalSAcorr_FA_glasser_hcpd) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_hcpd))) #set column names; 100 age x 10,000 posterior draw correlations df

deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws

cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")

ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
  geom_ribbon(fill = c(alpha("#323280", 0.4))) +
  geom_line(linewidth = .2, color = "black") +
  labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
  theme_classic() + 
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
  #scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
  theme(
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.line = element_line(linewidth = .2), 
   axis.ticks = element_line(linewidth = .2)) 

#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]

#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
  df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
  df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes,  by = c("label", "tract", "orig_parcelname"))
  
  #S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - development correlation
  anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")

anatomical.ps <- c(0.0001762558, 1.604827e-10, 0.000236177)
p.adjust(anatomical.ps, method = c("fdr"))

#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4859962)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3097694)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1192805)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = -0.09059357)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "D-V", "M-L"))

ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) + 
  geom_bar(stat='identity') + 
  labs(x="") +
  labs(y="\nAge of Maturation Correlation") +
  theme_classic()+
  scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_pnc.pdf", device = "pdf", dpi = 500, width = 1.4, height = 1)

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")

anatomical.ps <- c(0.3912361, 7.947443e-11, 4.340955e-06)
p.adjust(anatomical.ps, method = c("fdr"))

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)

#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)

#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full

library(ciftiTools)

#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels

#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version 
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released 
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)

#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")

#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

detach("package:ciftiTools", unload=TRUE)

EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)

df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = amplitude.agedecline, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(7.9, 18), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.75, height = 5.4) 

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c( "#581768", "#943b90", "#EA8B57", "#fcd16d"), limits = c(0.005, 0.018), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.44, height = 5.13) 

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(-0.04, -0.0345), oob = squish, na.value = "white") 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.13, height = 5.08) 

plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)

EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)

myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)

amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

plasticity.ps <- c(EI.thalamus.pnc, myelin.thalamus.pnc, amplitude.thalamus.pnc)
plasticity.ps

p.adjust(plasticity.ps, method = c("fdr"))

plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)

EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)

myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)

amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") 

ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

plasticity.ps <- c(EI.thalamus.hcpd, myelin.thalamus.hcpd, amplitude.thalamus.hcpd)
plasticity.ps

p.adjust(plasticity.ps, method = c("fdr"))